DEanalysis_kfish_osmotic

Introduction

A mixed model with random effects was chosen for this multifactor experiment, and analyzed using the limma package in R. This package implements a linear modeling approach and empirical Bayes statistics. The limma package with the voom method estimates the mean-variance relationship of the log-counts, generates a precision weight for each observation and enters these into the limma empirical Bayes analysis pipeline.

In this model, clade (levels = Clade1, Clade2, Clade3), physiology (levels = Marine, Freshwater), and experimental condition (levels = 15ppt, 0.2ppt) are fixed effects while species (levels = 14) is considered a random effect.

The raw counts file, generated with NumReads from the salmon (version 0.12.0) quantification tool and summarized with the tximport Bioconductor package (version 1.10.1) in R, can be downloaded from an osf repository with this link, then imported into the R framework.

Samples from species with low numbers of replicates were dropped from the raw counts table (F. zebrinus, F. nottii, F. sciadicus). The raw counts table has the following dimensions (genes x samples).

#dim(design)
#[1]  5 81
dim(counts)

[1] 30466    81

Sample Design Matrix

A matrix was generated using the following model with fixed effects:

~0 + physiology*condition*clade

The random effect of species will be taken into account later.

Filtering and Normalization

Genes with low expression across samples were dropped from the analysis using a conservative approach. The function filterByExpr was used on the raw counts matrix. For each condition_physiology group (regardless of species), each sample must have a minium count of 10, and a group minium total count of 100. This reduced the counts table to the following dimensions (genes x samples):

gene.names<-rownames(counts)
counts<-as.matrix(as.data.frame(sapply(counts, as.numeric)))
rownames(counts)<-gene.names
#class(counts)

keep<-filterByExpr(counts,design = design,group=condition_physiology,min.count = 10, min.total.count = 100)
counts.filt <- counts[keep,]
dim(counts.filt)
#write.csv(counts.filt,"../../../21k_counts_filt_30April2019.csv")

Genes of Interest

After filtration, we check the counts matrix for the presence of several genes of interest. These genes have demonstrated responsiveness to salinity change from past studies.

gene	Funhe	Ensembl
zymogen granule membrane protein 16	Funhe2EKm029929	ENSFHEP00000007220.1
zymogen granule membrane protein 16	Funhe2EKm029931	ENSFHEP00000025841
solute carrier family 12 member 3-like (removed)	Funhe2EKm006896	ENSFHEP00000009214
chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed)	Funhe2EKm024148	ENSFHEP00000019510
ATP-sensitive inward rectifier potassium channel 1	Funhe2EKm001965	ENSFHEP00000015383
inward rectifier potassium channel 2	Funhe2EKm023780	ENSFHEP00000009753
aquaporin-3		ENSFHEP00000006725
cftr	Funhe2EKm024501	ENSFHEP00000008393
polyamine-modulated factor 1-like	Funhe2EKm031049	ENSFHEP00000013324
sodium/potassium/calcium exchanger 2	Funhe2EKm025497	ENSFHEP00000034177
septin-2B isoform X2		ENSFHEP00000015765
CLOCK-interacting pacemaker-like	Funhe2EKm026846	ENSFHEP00000017303
vasopressin V2 receptor-like	Funhe2EKm026721	ENSFHEP00000000036
sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1	Funhe2EKm001174	ENSFHEP00000031108
septin-2	Funhe2EKm012182	ENSFHEP00000016853
otopetrin-2	Funhe2EKm035427	ENSFHEP00000026411
claudin 8	Funhe2EKm037718	ENSFHEP00000006282
claudin 4	Funhe2EKm007149	ENSFHEP00000003908

If the Ensembl ID displays below, the gene is present in the whole data set and has not been filtered.

# ---------------------------
# Andrew Whitehead's genes of interest:
# ---------------------------

# zymogen granule membrane protein 16
# Funhe2EKm029929
# ENSFHEP00000007220.1
countsfilt <- counts.filt
countsfilt$row <- rownames(countsfilt)
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000007220.1"]
goi
# zymogen granule membrane protein 16
# Funhe2EKm029931
# ENSFHEP00000025841
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000025841"]
goi
# solute carrier family 12 member 3-like (removed) 
# Funhe2EKm006896
# ENSFHEP00000009214
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009214"]
goi
# chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed)
# Funhe2EKm024148
# ENSFHEP00000019510
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000019510"]
goi
# ATP-sensitive inward rectifier potassium channel 1 
# Funhe2EKm001965
# ENSFHEP00000015383
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015383"]
goi
# inward rectifier potassium channel 2
# Funhe2EKm023780
# ENSFHEP00000009753
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009753"]

# --------------------------------
# other salinity genes of interest
# --------------------------------
# ============================================
# aquaporin-3
# ENSFHEP00000006725
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006725"]
goi
# cftr
# Funhe2EKm024501
# ENSFHEP00000008393
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000008393"]
goi
# polyamine-modulated factor 1-like
# Funhe2EKm031049
# ENSFHEP00000013324
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000013324"]
goi
# sodium/potassium/calcium exchanger 2
# ENSFHEP00000034177
# Funhe2EKm025497
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000034177"]
goi
# septin-2B isoform X2
# ENSFHEP00000015765
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015765"]
goi
# CLOCK-interacting pacemaker-like
# ENSFHEP00000017303
# Funhe2EKm026846
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000017303"]
goi
# vasopressin V2 receptor-like
# Funhe2EKm026721
# ENSFHEP00000000036
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000000036"]
goi
# sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1
# ENSFHEP00000031108
# Funhe2EKm001174
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000031108"]
goi
# septin-2
# Funhe2EKm012182
# ENSFHEP00000016853
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000016853"]
goi
# otopetrin-2
# Funhe2EKm035427
# ENSFHEP00000026411
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000026411"]
goi
# claudin 8
# Funhe2EKm037718
# ENSFHEP00000006282
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006282"]
goi
# claudin 4
# ENSFHEP00000003908
# Funhe2EKm007149
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000003908"]
goi
all_goi<-c("ENSFHEP00000007220.1","ENSFHEP00000025841","ENSFHEP00000019510",
           "ENSFHEP00000015383","ENSFHEP00000009753","ENSFHEP00000006725","ENSFHEP00000008393",
           "ENSFHEP00000013324","ENSFHEP00000001609","ENSFHEP00000013324","ENSFHEP00000034177",
           "ENSFHEP00000015765","ENSFHEP00000017303","ENSFHEP00000000036","ENSFHEP00000031108",
           "ENSFHEP00000016853","ENSFHEP00000003908")

Exploratory Plots

Log counts before normalization:

# log counts before DE
boxplot(log(counts.filt+1), las = 2, main = "")

Log counts after cpm normalization

# ---------------

# DE analysis

# ---------------

genes = DGEList(count = counts.filt, group = condition_physiology_clade)
genes = calcNormFactors( genes )

# write normalized counts
dir <- "~/Documents/UCDavis/Whitehead/"
tmp <- as.data.frame(cpm(genes))
tmp$Ensembl <- rownames(tmp)
tmp <- dplyr::select(tmp, Ensembl, everything())
write.csv(tmp, file = file.path(dir, "normalized_counts.csv"), quote = F, row.names = F)

vobj = voom( genes, design, plot=TRUE)
lcpm <- cpm(genes$counts, log = TRUE)

# log counts after DE

boxplot(lcpm, las = 2, main = "")
plot(colSums(t(lcpm)))

Voom weights:

vwts <- voomWithQualityWeights(genes, design=design, normalization="quantile", plot=TRUE)

The random effects of species are taken into account with the duplicateCorrelation function, which estimates the correlation between species. This reflects the between-species variability. The higher the correlation (0-1.0), the higher the variability between species.

corfit <- duplicateCorrelation(vobj,design,block=ExpDesign$species)

corfit$consensus
#[1] 0.751381

Individuals clustered by overall expression

Individuals clustered by Top 100 genes

Individuals clustered by top 50 gene expression

PCA for overall expression

Contrasts

After running the lmFit function, which fits the linear model for each gene in the matrix and takes the random effects correlation into account, the resulting linear model fit is then used to compute the estimated coefficients and standard errors for a given set of contrasts:

fitRan <- lmFit(vobj,design,block=ExpDesign$species,correlation=corfit$consensus)
#colnames(coef(fitRan))
y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|clade*physiology ), reverse = T)
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$clade, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts
contrasts

[1] "15_ppt_v_0.2_ppt_Clade1_FW" "15_ppt_v_0.2_ppt_Clade2_FW" "15_ppt_v_0.2_ppt_Clade3_FW"
[4] "15_ppt_v_0.2_ppt_Clade1_M"  "15_ppt_v_0.2_ppt_Clade2_M"  "15_ppt_v_0.2_ppt_Clade3_M" 

tables <- lapply(contrasts, function(contr){
    print(contr)
    cm <- contrast.matrix[contr,]
    ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
    cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
    tmp <- contrasts.fit(fitRan, contrasts = cm)
    tmp <- eBayes(tmp)
    tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
    tmp3 <- tmp2
    tmp3$row <- rownames(tmp3)
    tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
    tmp3 <- tmp3[order(tmp3$adj.P.Val),]
    filename <- paste0(contr, ".csv")
    #write.csv(tmp2, file = file.path(dir, filename),quote = F)
    tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
    header1 <- 6
    names(header1) <- paste0("Top 20 genes for ", contr)
    header2 <- 6
    names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
    header3 <- 6
    names(header3) <- paste0("Output file is ", filename)
    tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
    tab <- tab %>% kable_styling()
    return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}
)

[1] "15_ppt_v_0.2_ppt_Clade1_FW"

row names of contrasts don't match col names of coefficients

[1] "15_ppt_v_0.2_ppt_Clade2_FW"

row names of contrasts don't match col names of coefficients

[1] "15_ppt_v_0.2_ppt_Clade3_FW"

row names of contrasts don't match col names of coefficients

[1] "15_ppt_v_0.2_ppt_Clade1_M"

row names of contrasts don't match col names of coefficients

[1] "15_ppt_v_0.2_ppt_Clade2_M"

row names of contrasts don't match col names of coefficients

[1] "15_ppt_v_0.2_ppt_Clade3_M"

row names of contrasts don't match col names of coefficients

Three-way contrasts

The number of genes significant for the three-way interaction of condition:physiology:clade:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")

Two-way contrasts

These genes demonstrate a conserved response to experimental condition across M or FW physiologies, regardless of clade.

y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|physiology ), reverse = T)
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts

contrasts

tables <- lapply(contrasts, function(contr){
  #print(contr)
  cm <- contrast.matrix[contr,]
  ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
  cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
  tmp <- contrasts.fit(fitRan, contrasts = cm)
  tmp <- eBayes(tmp)
  tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
  #tmp3 <- tmp2
  #tmp3$row <- rownames(tmp3)
  #tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
  #tmp3 <- tmp3[order(tmp3$adj.P.Val),]
  filename <- paste0(contr, ".csv")
  #write.csv(tmp2, file = file.path(dir, filename),quote = F)
  tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
  header1 <- 6
  names(header1) <- paste0("Top 20 genes for ", contr)
  header2 <- 6
  names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
  header3 <- 6
  names(header3) <- paste0("Output file is ", filename)
  tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
  tab <- tab %>% kable_styling()
  return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}

)

The number of genes significant for the two-way interaction of condition:physiology, independent of clade:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")

Condition - main effect

These genes demonstrate a main effect of condition (15 ppt vs. 0.2 ppt), regardless of clade or physiology.

y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition), reverse = T)
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
rownames(contrast.matrix) <- contrasts

contrasts

tables <- lapply(contrasts, function(contr){
  #print(contr)
  cm <- contrast.matrix[contr,]
  ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
  cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
  tmp <- contrasts.fit(fitRan, contrasts = cm)
  tmp <- eBayes(tmp)
  tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
  #tmp3 <- tmp2
  #tmp3$row <- rownames(tmp3)
  #tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
  #tmp3 <- tmp3[order(tmp3$adj.P.Val),]
  filename <- paste0(contr, ".csv")
  #write.csv(tmp2, file = file.path(dir, filename),quote = F)
  tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
  header1 <- 6
  names(header1) <- paste0("Top 20 genes for ", contr)
  header2 <- 6
  names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
  header3 <- 6
  names(header3) <- paste0("Output file is ", filename)
  tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
  tab <- tab %>% kable_styling()
  return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}

)

The number of genes significant for the main effect condition:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")

Heatmaps

Marine-Clade 1 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Freshwater-Clade 1 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Marine-Clade 2 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Freshwater-Clade 2 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Marine-Clade 3 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Freshwater-Clade 3 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Marine physiology (two-way) conserved response

rld <- log2(mean_norm_counts_ordered_M_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Freshwater physiology (two-way) conserved response

rld <- log2(mean_norm_counts_ordered_FW_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

Condition response (main effect)

rld <- log2(mean_norm_counts_ordered_condition_sig+1)
geneDists <- dist(mean_norm_counts_ordered_condition_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

sessionInfo()